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In a granular solid, mechanical equilibrium requires a delicate balance of forces at the disordered 
grain scale. To understand how macroscopic rigidity can emerge in this amorphous solid, it is crucial 
that we understand how Newton's laws pass from the disordered grain scale to the laboratory scale. 
In this work, we introduce an exact discrete calculus, in which Newton's laws appear as differential 
relations at the scale of a single grain. Using this calculus, we introduce gauge variables which 
describe identically force- and torque-balanced configurations. In a first, intrinsic formulation, we 
use the topology of the contact network, but not its geometry. In a second, extrinsic formulation, we 
introduce geometry with the Delaunay triangulation. These formulations show, with exact methods, 
how topology and geometry in a disordered medium are related by constraints. In particular, we 
derive Airy's expression for a divergence-free, symmetric stress tensor in two and three dimensions. 



INTRODUCTION 

Poised between fluids and solids, granular media pose 
outstanding challenges in fundamental physics [![. At 
finite pressure, a cohesionless, frictional granular medium 
resists bulk and shear deformation, and is therefore a 
solid. However, application of a sufficiently large shear 
stress will cause the medium to deform indefinitely. It 
is largely unknown how grain-scale quantities vary with 
the applied stress, and therefore we cannot yet predict 
the yield stress from first principles. 

At the grain scale, the mechanical stability of the solid 
phase requires a delicate balance of forces and geometry. 
For nearly rigid grains, like sand, a good approximation 
is that intergranular forces are transmitted at discrete 
contacts. The positions of these contacts, and the forces 
across them, must be arranged so that the net force and 
net torque on each grain vanish. To understand how 
a solid granular medium responds to a macroscopically 
applied stress, we need to understand how these "micro- 
scopic" constraints pass from the disordered grain scale 
to the laboratory scale. 

To provide insight on this process, in this note we show 
how forces and torques can be derived from discrete po- 
tentials which identically satisfy these constraints. These 
potentials have continuous analogs which can be exactly 
related to grain-scale quantities by a discrete calculus, 
which we develop. In particular, we consider a static, 
frictional packing f2 in dimensions d = 2 and d = 3, in 
the absence of body forces. Each of the N soft, repulsive, 
identical disks (d = 2) or spheres (d = 3) g is subject to 
force and torque balance: 

= E ^> = E ( rC - r3 ) x n- 

Here f£ is the contact force exerted on grain g at contact 
c, r c is the position of c, and r 9 is position of the center 



of g. Each contact force must also satisfy a repulsive 
condition (r c — r 9 ) ■ /° < 0. The macroscopic object of 
interest is the stress tensor (ij 

gGGceCs 

where V is the volume of the packing. In what follows we 
will decompose this into contributions from each grain g, 
a 9 = 1/V 9 J2 ce cs( rC ~ r9 )fgi where V 9 is a volume 
associated to a grain, discussed below. 

If a continuum description exists, then mechanical 
equilibrium requires that the stress tensor satisfies Q 

= V • <r, <r = & T . (3) 

This is the macroscopic analog of ([T]) . In two dimensions 
(2D), in the absence of body forces, it implies that cr can 
be written as Q 

6- = V x V x ip, (4) 

where ip is a scalar, known as the Airy stress function. 
In component form, 

& = f dyyip -d xy ip\ 

\-dxy1p dxxlp J ' 

In three dimensions (3D), (j4|) also holds, with ip replaced 
by a tensor xp, known as the Beltrami stress tensor. A <r 
of this form is identically divergence-free and symmetric, 
so ip automatically yields stress configurations in force 
and torque balance. The extension of (j4|) to granular 
materials was initiated by the seminal works of Satakc 
[H and Ball and Blumenfeld jf|. These authors found 
"loop forces" p such that <r = V x p, as well as par- 
tial analogs of ip. However, they were not able to derive 
p from a discrete version of ip which gives identically 
torque-balanced configurations, so that the analogy with 
dU is incomplete. 
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In this paper we extend their results by deriving an ex- 
act, discrete analog of Q for frictional sphere packings in 
two and three dimensions. We first present an intrinsic 
formulation based solely on the topology of the contact 
network. We then present a complementary extrinsic for- 
mulation, by building the contact network into a triangu- 
lation of space. The added structure of the triangulation 
depends heavily on the geometry of the contact network. 

The new variables p and ip define changes of coordi- 
nates in phase space. In this paper, we consider the limit 
of large grain rigidity in which the deformation of grains 
can be ignored. This allows us to fix the grain positions 
and consider the Force Network Ensemble of force con- 
figurations on a fixed packing geometry jjjj]. 

The geometry consists of Nrg real (force-bearing) 
grains as well as Nyo virtual grains, or rattlers, which are 
trapped in the packing, but do not contribute to mechan- 
ical stability. The real grains touch at Nrc force-bearing 
contacts § , so phase space is spanned by the Nrc con- 
tact forces and has dimension <INrc- Force and torque 
balance restrict force configurations to a subset of phase 
space, of dimension 



M d = alN RC 



d(d + l) 



N 



RG 



(5) 



Writing Nrc = zNrq/2, wc define the mean contact 
number z. When Md = 0, the forces are uniquely deter- 
mined by the grain positions; this is known as isostaticity 

0, [l(| and occurs when z = d + 1. More generally we 
may also consider hyperstatic packings with z > d + 1 
and Md > 0. 

Since Md depends only on the geometry, it is a topo- 
logical invariant under any coordinate change in force 
space. The intrinsic formulation uses the contact net- 
work, which assigns a vertex to each grain center and a 
link to each contact. Mechanical stability requires that 
real grains form closed loops I. This property implies a 
combinatorial identity between the number of indepen- 
dent loops Nl, real grains Nrg, and real contacts Nrc, 
which we can use to rewrite Md- To see this identity, 
we consider an arbitrary grain g and inductively build 
the entire packing, one loop at a time. Initially, we have 
one grain, no contacts, and no loops. We choose an arbi- 
trary contact belonging to g and trace out a shortest 
loop back to g, containing, say, k new grains. In doing 
so we have added k grains, k + 1 contacts, and 1 loop, 
so Nrg — Nrc + Nl is preserved at its original value, 

1. Continuing in this way, always beginning loops on ex- 
isting grains, we eventually trace out the entire interior 
contact network, so that 



N, 
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RG 
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This is a version of Euler's formula [12j for the contact 
network. We call the set of independent loops chosen by 
this process a net. In 2D, the net is unique, but in higher 
dimensions, many nets are possible (l3| . 



The relation (|6]) allows us to write M d in two distin- 
guished ways, by eliminating either Nrc, or Nrg- The 
first, 



M d = dN L - ^L-^Nr G 



(7) 



indicates that phase space is spanned by a vector field 
defined on the loops, providing dN^ degrees of freedom, 
subject to torque balance, providing \d{d — 1)Nrg con- 
straints, and a single extra vector constraint. To exhibit 
this formulation explicitly, we assign to each loop a fixed 
orientation, and a pseudovector loop force p e . Setting 
p e g = +p e if the oriented loop £ enters the grain g, and 
p g = p e if it exits g, wc now write each contact force 
as 



/ S C =E^ 



(8) 



where L c is the set of loops adjacent to the contact c. 
This definition ensures that Newton's 3rd law is satis- 
fied. Moreover, since each loop going through g enters at 
precisely one contact, and exits at precisely one contact, 
when we sum the contact forces incident on a grain, each 
loop force appears in equal and opposite pairs, so that 
force balance is satisfied identically. 

In 2D, each contact is adjacent to two loops, and the 
loops can all be taken to have anticlockwise orientation. 
This makes |(5J) a simple difference of loop forces, so that 
adding a constant to each p e leaves invariant the contact 
forces. In this case, there is a vector gauge freedom and 
so a vector constraint is needed to fix a gauge. However, 
in higher dimensions, the net is not unique, and there is 
no natural orientation of the loops. It is not clear in this 
case whether a gauge freedom exists, or whether the extra 
constraint is a net-dependent consistency requirement. 

We may also write 



M d 



d{d + r 



Nr 



d(d-l) 



N 



RC 



d(d+l) 



(9) 



in which the number of grains no longer appears. This 
formulation is achieved by supplementing p e with a field 
ip e , a pscudoscalar in 2D and a pseudovector in 3D, de- 
fined so that 



(10) 



where r e is the position of i, discussed in the sequel. As 
with p^, (p e — +(p e if the oriented loop i enters the grain 
g, and ip g = —ip e if it exits g. Summing the torques 
applied to a grain, it is easily seen that torque balance is 
satisfied identically. 

To ensure that torques result from tangential contact 
forces, the torques inferred from p l must equal those de- 
termined from cp e , which requires 
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-r c ) 



(11) 
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These are the contact constraints appearing in ((9j)- As 
with the p formulation, the residual constraint may de- 
pend on the net. In two dimensions © and (JTUJ) reduce 
to Satake's formulation 

The intrinsic formulation can be put into a differential 
form by defining topological divergence and curl opera- 
tors. The former, defined for vector fields on the con- 
tacts, is (d*f) 9 = J^cecs fg> while the latter, defined 
for pseudovector and pscudoscalar fields on the loops, is 
{d*p) c g = J2eeL c Pg- These operators satisfy the identity 
d*d* = 0. Dropping subscripts and superscripts, we can 
write Newton's laws as d*f = and d*(r x /) = 0, ([£} as 
/ = d*p, and (dHJ as r x / = d*(tp + r x p). 

Continuing in this fashion, one can define a complete 
topological calculus, with analogs of Stokes Theorem and 
Poincare's Lemma. There is a difficulty, however; vec- 
tor fields, such as /, have no direct continuum meaning, 
since they change sign with the orientation of the con- 
tact. They are, properly, 1-forms. This difficulty may 
be overcome by resolving / along a vector field; this is 
precisely what is accomplished by the stress tensor, ([2]). 
This introduces geometry into the formulation. In 2D 
the loops tile the plane, so their geometry is simple, but 
in 3D, they may connect in a complex way. Moreover, in 
3D nets are not uniquely defined, so given a contact c, 
there is no way of knowing how many loops are adjacent 
to c without some knowledge of the whole net. 

Rather than proceeding with this admixture of topol- 
ogy and geometry, we therefore present a complementary 
extrinsic formulation based on a triangulation of space. 
By using a triangulation, the topology is locally regular 
and simply related to the geometry. Around any grain, 
we can determine the topology and geometry of the trian- 
gulation knowing only the positions of the grain's neigh- 
bours. In particular, we use the Delaunay triangulation 
and its dual, the Voronoi tessellation [l4[, shown in Fig- 
ure [1] 

The Delaunay triangulation is formed by adding to the 
contact network links between adjacent grains which do 
not touch, so that space is partitioned into simplices: 
triangles in 2-space and tetrahedra in 3-space [l5| . Phys- 
ically, the new "virtual" contacts correspond to contacts 
that could be created under a small deformation flil ]. 
The Voronoi tesselation is the dual of the Delaunay tri- 
angulation; each Voronoi cell is centered on a grain and 
contains the points nearer to that grain than any other. 
We define the volume associated to a grain, V 9 , as the 
volume of its Voronoi cell. 

Our general approach to derive an analog of (J4j) is to 
define an exact discrete calculus 17J on these graphs in 
which (JTJ appear as differential relations. By Poincare's 
Lemma, V-F = 0=>F = VxG, these suggest natural 
forms for p and ip. Routine calculations then show that 
(TTJ) holds identically in the new variables. 



Notation 

In this paper we deal with scalars, vectors, tensors, and 
their pseudo counterparts, defined on grains g, contacts c, 
triangles t, and tetrahedra v. When unambiguous, these 
are indicated by a single superscript, but it is often nec- 
essary to add subscripts to indicate orientation. These 
are explained as they appear. 

The sets of all real grains, virtual grains, real contacts, 
virtual contacts, triangles, and tetrahedra in the packing 
are denoted by RG, VG, RC , VC, T, and V, respec- 
tively. Together, the real and virtual grains form the set 
of grains G = RGUVG, and similarly, the real and virtual 
contacts together form the set of contacts C = RCUVC. 
We also consider local sets, for example C 9 denotes the 
set of contacts surrounding a grain g, and likewise for 
other neighboring quantities. The sets of boundary con- 
tacts, triangles, and tetrahedra are denoted BC, BT, and 
BV, respectively. 

Tensors are denoted with a hat. The identity tensor 
is written 5. Outer (dyadic) products are simply de- 
noted with a space. For vectors and 2-tensors, opera- 
tors always act on adjacent indices, in the same order 
as they are normally written, and all contractions are 
indicated with a dot, the number of dots indicating the 
number of contracted indices; e.g., A : X7(B (V ■ C)) = 
Aijdi(Bj k (diCi m )). We only need one 3-tensor, 
the 3D Lcvi-Civita symbol appearing in the 3D cross 
product. We take the free index to be the first, e.g., 
(A x B)ij m = J2k,iAik£jidBi m . To conform to conven- 
tion we make an exception to the above rules for the 
tensor curl in 3D, letting (V x B)^ = J2 k -,;./'>;. 

In 2D the cross product is, e.g., (A x B)u = 
y~]. k AijSjkBki, with f the 2D Levi-Civita symbol, £12 = 
— e 2 i = 1, £11 = £22 = 0. The 2D curl is (V x u)ij = 
J2k £ ikdkUj. In both 2D and 3D, the condition that a 
2-tensor A be symmetric can be written 



A = A 



i : A = 0. 



(12) 



Contours are always oriented anticlockwise, and we 
use +/— to indicate geometric elements anticlock- 
wise/clockwise from given elements, around grains, tri- 
angles, etc. For example, in 2D, c + = c + (g,i) is the 
contact anticlockwise from triangle t when looking from 
grain g. 



2D 



In the plane, contacts correspond to edges of Voronoi 
cells; we assign vectors s c g along these edges, circulating 
anticlockwise around grains (Figure [TJ . We also assign 
vectors i g pointing from the center of g to the center of 
its neighbour at c. These allow a natural definition for 
the area of a contact, A c = | £ c x s c . For a tensor field 



FIG. 1. Delaunay triangulation (left) and Voronoi tesselation (right) in 2D. Real contacts are indicated by small dots, with their 
associated contact vectors solid. Virtual contacts are shown as unfilled circles, with their associated contact vectors dashed. 
Although in reality soft spheres deform on contact, for simplicity real contacts are shown with overlaps. 



defined on the contacts, we define its divergence by anal- 
ogy with a tensor version of Gauss' Theorem, defining 
the resulting line integral with the s vectors: 



/ dAv -a = A 9 (v ■ a) 9 = - s c g xa c = -<£ dr 



We emphasize that the underlying fields and definitions 
are discrete. Because the discrete fields satisfy a similar 
algebra as in the continuum, the continuum notation is 
useful for intuition and calculations. 

Rewriting @ as a sum over contacts, we have A c a c = 
—£ c f c . From this definition we see that s c g x a c = 2f c , 
so (V • a) 9 = is equivalent to force balance. By in- 
spection of ©, and using (|12[) . we see that a 9 = (a 9 ) T 
is equivalent to torque balance, and hence we reproduce 
the expected continuum equations ([3]) in the discrete cal- 
culus, at the grain scale. 

We will define the pseudovector p on the triangles, 
which correspond to vertices of Voronoi cells; we assign 
vectors s* circulating anticlockwise around grains, con- 
necting the Voronoi vertices. Using a form of Green's 
Theorem, for a vector field defined on the triangles, we 
may define its curl as 

f dAV x p = A 9 (V x pf = - J2 4 P* = - I dr p. 



teTs 



(13) 

Noting that 2s* = if - £ c g , we see that <r 9 = (V x p) 9 
if 

/ 9 C = P 4 ~-P t+ , (14) 

where t + = t + (g,c) is the triangle anticlockwise from 
contact c. Summing the contact forces incident on a 
grain, all loop forces appear in equal and opposite pairs, 
so force balance is identically satisfied. 



By working with the Delaunay triangulation, we treat 
real and virtual contacts on the same footing. Therefore, 
equation (|14[) applies to both real and virtual contacts, 
and hence generically leads to virtual contact forces. To 
obtain physical force configurations, we must explicitly 
impose that these vanish: if c e V C 9 , = f g = p* — 

p t+ . Since virtual contacts only occur in the interior of 
loops, this shows that p* must be constant on loops; this 
also shows that we recover the formulation of Satake Q 
and Ball and Blumenfeld [f| if we sum over the virtual 
contacts. 

An immediate consequence of ff 9 = (Vx p) 9 is that the 
stress tensor for the packing can be written as a boundary 
sum [21 



dA a = Aa 



E 



dr p, 



nil 



where the s's connect boundary contacts in an anticlock- 
wise manner. 

We may check that is conserved under this coor- 
dinate change by counting triangles. Applying Euler's 
formula to a single loop £, Ny G — Ny C + = 1. 
Adding this relation over the whole packing, we find 
N VG - N vc + N T = N L , so that, using ©, M 2 = 
2N T - (2N VC + N RG - 2N VG + 2). This indicates that 
the loop forces are constrained by the virtual contact 
constraints and by torque balance, and that 2Ny G con- 
straints become redundant. This redundancy is a con- 
sequence of the definition of loop forces. Indeed, since 
loop forces automatically yield grains in force balance, 
for a virtual grain with z virtual contacts, it is sufficient 
to impose z — 1 virtual contact contraints to guarantee 
that all z virtual contact forces vanish. Finally, the loop 
forces have a gauge freedom p —> p + Ap, with Ap any 
constant, as is clear from their definition. 
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The torque balance constraint can be rewritten in 
terms of p by first defining 



dA Vp = A 9 (V • p) 



E4xp* 



Comparing this equation with (fT3| , and using & 9 = ( V x 
p) 9 , we see by inspection that <r 9 = (<r 9 ) T is equivalent 
to (V • p) 9 = 0. This motivates a search for ip c such that 
p* = V x ip c . We define 



dr ip. 



f dAV x ip = A 1 {V x ip) 1 = - ^ t t i) c = - I , 

cgct Jet 



It is not obvious from this definition, but as shown in 
Appendix A „ this yields force configurations that iden- 
tically satisfy torque balance. Since each triangle is sur- 
rounded by three contacts, and contacts interior to a loop 
are shared by two triangles, = 2Ny C + N RC . This 
yields the global counting relation 3Nt = 2Nc + Nbc, 
in which BC is the set of triangle edges that circulate 
around the boundary. We can therefore write Mi = 
(N c + N BC ) - (2N VC - 3A yG + 3 )- The stress fac- 
tion ip is defined on the real and virtual contacts, and 
constrained on the virtual contacts by the requirement 
of no virtual contact force. Because rattlers automati- 
cally satisfy force and torque balance, 3Nvg constraints 
are redundant. Finally, ip has a three dimensional gauge 
freedom ip — > ip + Aip, with (Aip) c = A-r c + B any linear 
function of position, which follows from the computation 

V x (AV>) C = £■ A. 

We have therefore succeeded in writing <r 9 = V x p t = 

V x V x ip c . It is possible to invert this transformation, 
up to gauge freedom. For a contour of adjacent triangles 
(to, t\, . . . , t n ), separated by contacts (ci, C2, . . . , c„), we 
find 



(15) 

where the forces are exerted from the left side of the 
contour to the right side. Similarly, 



n-l 



V> c - - ip c 



i=0 



dr x p. (16) 



Force and torque balance ensure that the line integrals 
in (fl5|) and (|T6|) . respectively, are path- independent. 

In the 2D discrete calculus, we therefore have the the- 
orems 



and 



V • <r c = 



V-p* 



& 9 = v x p l 



(17) 



(18) 



Equation (|T6|) allows us to relate ip c to ip in the intrin- 
sic formulation. We first allow ip to take two values on 
each real contact, one per adjacent loop ±£ We define a 



dry. p mean ip within a loop by z ip = J2 c eRC e wnere 



z is 

the number of contacts around the loop. Then (|16[) im- 
plies ipl = ip e + (r e — r c ) x p £ , where z l r l = ^ c£RC i r c . 
Comparing with (fTTj) we see that cp e = ip e and the consis- 
tency constraints in the intrinsic formulation are equiva- 
lent to ip C j ri = ip°_i, a type of Newton's 3rd law. 




FIG. 2. Subset of Delaunay triangulation and Voronoi tesse- 
lation in 3D. Real contacts are indicated by small dots, with 
their associated contact vectors solid. Virtual contacts are 
shown as unfilled circles, with their associated contact vec- 
tors dashed. 



3D 



In three dimensions, we can proceed analogously. The 
Voronoi tesselation assigns to each grain a convex poly- 
hedron (Figure^]). Each face corresponds to a contact, 
each edge corresponds to a Delaunay triangle, and each 
vertex corresponds to a Delaunay tetrahedron. As before, 
we can define a discrete divergence with Gauss' Theorem: 

f dv v ■ a = v 9 (v • a) 9 = J2 A % &c = [ dS a > 

where £ = £/\£\ and the natural area for a contact, A c . is 
the area of the contact face. Rewriting ([2]) as a sum over 
contacts, we have V c a c = —£ c f c . The natural choice of 
V c is the volume of the bipyramid whose vertices are the 
neighbouring grain centers and whose base is the contact 
face. Using V c = ±A C \£ C \, we see that A c £ c g ■ <r c = -U c g ■ 
£ c f c = — 3/g, so again (V ■ a) 9 = is equivalent to 
force balance. Again a 9 = (a 9 ) T is equivalent to torque 
balance, and hence we reproduce the expected continuum 
equations @ in the discrete calculus, at the grain scale. 
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We now seek to define p so that <r = V x p. By Gauss' 
Theorem 

f dV Vxp = V 9 (V x p) 9 = J2 dS*x(p*) T = / dSxp 

J a ^*-<n„ J da 



The natural way to define dS g is to decompose the con- 
tact faces on either side of t into triangles of area A l c ± and 
t pc- = J2 cect Alfcg. It is easily 



set dS l g 



A* if +A t __£ c „ 



seen from Figurc[2]that A^i^ = — r c ) x s gc , where r* 
is the intersection of the triangle t with the corresponding 
edge of the Voronoi cell of g, and s gc circulates anticlock- 
wise around contacts. We can now rewrite (V x p) 9 as 

l"(Vxpf = 4EE ((r c - r f ) x si) x (p'f 



4EE r-r 9 )x S ; c )x(^ 
-IE E(^x4c)x(pT, 



where we have used the identity X^cec* s gc = to ex- 
change r* with r 9 . Using ([2]), we now see that cr 9 = 

v x p* if ^/ ff c = ±£ teTC x 4c) x OT- This is 

equivalent to 



"ffc — 2 



i6T<= 



(<)c 



pdr, (19) 



and = X)teT c s gcP* ' ^o- The latter equation is satisfied 
identically if we let p* = — p gc s gc - This allows (fT9|) 
to be rewritten f g = X)teT<= /^gc which shows explicitly 
that the contact forces depend only on a vectorial quan- 
tity, albeit one with an intrinsic orientation. It also shows 
reduction to the form of the intrinsic formulation. 

Summing the contact forces incident on a grain, all 
loop forces appear in equal and opposite pairs, so force 
balance is identically satisfied. As in 2D, we must ex- 
plicitly require that there be no virtual contact forces. A 



new feature in 3D is that p gc has a nontrivial gauge free- 
dom p gc -> p gc + (ApY gc with (Ap)* c = B v+ - B" , for 
any vector field B v defined on the Delaunay tetrahedra, 
which are dual to Voronoi vertices. 

Again, ^ 9 = (Vx p) 9 implies that the stress tensor for 
the packing can be written as a boundary sum Q 

f dV & = V&= A c l c x (p') T = <f dSx p T , 
Jn ceBC Jdn 

where the Z's are oriented outwards. 

We now seek -0 C such that p* = V x %p c . Using Stokes' 
Theorem, we set 

J(Vxip)-dS = A 1 (V x = ^ C - £C = J $' dr > 



where s = s/\s\ and the contour is oriented anticlock- 
wise around s*. A natural choice to guarantee p* = 
V x ip c is to set ip c = j£ c l c ip c ; then we have A t p t gc = 
— ^1^*1^2 c(£C t £g t ^P c , in close analogy with the 2D case. 
The computation that this leads to torque-balanced force 
configurations is almost identical to the 2D case, and for 
the same reasons, works only for identical spheres, tjj has 
a complex gauge freedom, discussed in Appendix B „ 



DISCUSSION AND EXTENSIONS 
Physical Interpretation 

As gauge variables, p and tp do not admit a unique 
physical interpretation. However, in 2D, the inversion 
formulae suggest a macroscopic interpretation of poten- 
tial differences. In particular, equation (|15[) indicates 
that a difference of loop forces gives the flux of stress 
across any contour connecting two triangles. Similarly, 
equation (|T5)) indicates that a difference of stress func- 
tion gives the flux of loop force across any contour con- 
necting two contacts. A linear change in stress function 
gives a constant flux of loop force, and hence no stress. 
Therefore, stresses correspond to curvature of the stress 
function. We prove in Appendix C that in 2D, in re- 
gions where macroscopic normal stresses are repulsive, 
the stress function is convex. 

By considering how a single contact force is written 
in terms of ip, it is also possible to obtain a microscopic 
interpretation. In 2D, an arbitrary contact force f g can 
be decomposed as 



E 



1 



A t(c<) 



(20) 



where C c is the set of contacts surrounding the loops 
adjacent to c, and the £ c are oriented in the same direc- 
tion as £ g . Considering rigid grains so that the geome- 
try is fixed, we see from this expression that an increase 
in i\) c directly increases the normal force at c. It also 
propagates to neighboring contacts, adding to their con- 
tact forces a component in the direction of £ g . If the 
neighboring loops contain virtual contacts, the distur- 
bance further propagates to the surrounding contacts in 
such a way that force and torque balance are preserved. 
We can therefore associate a change in ip c with a change 
in the normal force at c, as well as the necessary changes 
in normal and tangential forces at neighbours of c so as 
to preserve force- and torque-balance. 

In 2D, one can also obtain some geometrical intuition 
for these quantities by plotting the loop forces as points 
in force space, and drawing lines connecting the loop 
forces of adjacent loops, which are the contact forces. 
The Maxwell-Cremona reciprocal diagram thus obtained 
[l8[ is a tiling of polygons, each tile corresponding to a 
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FIG. 3. Subset of a Maxwell-Cremona reciprocal diagram. 
Each polygon corresponds to a single grain. Also shown is 
the subdivision into triangles for one particular grain, and a 
triangle corresponding to an external contact. 



single grain (Figure [3]) . The total area of the tiling de- 
pends only on the boundary loop forces, and hence only 
on the boundary contact forces. It can be estimated for 
large packings, as follows. 

We define nominal tile centers for each grain, p 9 , and 
divide each tile into triangles, one for each contact. The 
signed area of each triangle is a c g = \{p l — p 9 ) x (p e — 
p 9 ). For interior contacts, adding the contributions from 



each adjacent grain, we have 



a = a g+ 



(21) 



cGRC 

-i(e-d--e) : A& 1 
Adet(<x), 



The total area of the tiling is 



p 9 , the interior term is exactly 

£ a* =%(£■&■£): £ f c g £ c g (23) 
cei?c 

(24) 
(25) 

so that A = Adet(a) + 0(y~N). For periodic packings, 
the boundary term vanishes identically fl8j | . 

The above choice of p 9 does not guarantee that each 
a c is positive, which might be desirable in applications 
[lij . If we choose, instead, p 9 = det(<r)/P £ ■ r 9 , 
where P = 1/2 tr(tr) is the pressure, then we find 
a c = — det(<x)/(2P) f g _ ■ £ g _. In this representation, 
repulsion of the contact forces is equivalent to positivity 
of the areas. Moreover, X^ce-RC flC = Adctft) as before, 
so the boundary term must again be negligible for large 
enough packings. 

Since the tiling area depends only on the bound- 
ary forces, it is invariant under force-balance preserving 
changes of interior forces. In our formulation, localized 
force rearrangements can be generated by changing the 
stress function ip c at a particular contact, which changes 
the loop forces at adjacent loops. This construction can 
therefore be used to explore the Force Network Ensem- 
ble, and serves the same purpose as the "wheel moves" 
of Tighe and collaborators 3, HI ■ 

Attempts to extend the reciprocal diag ram construc- 
tion to 3D have not yet been successful [1S|] . In fact, since 
the loop forces are proper pseudovectors in 3D, they can- 
not simply be plotted as points in space and therefore 
cannot define the vertices of any polyhedra. Hence, if an 
analog of the reciprocal diagram exists, it is not obviously 
related to the loop forces defined in this work. 



^=E^E ( 22 ) 

ceRC ceEC 

where the second term sums over the Nec external con- 
tacts at the boundary, to correct for overcounting in the 
first sum. It is clear from Figure[3]that A is independent 
of the choice of p 9 , but the relative size of the terms in 
(|2"2")l is sensitive to this choice. For example, the trivial 
choice p 9 = makes all a c = 0, and the entire contribu- 
tion comes from the boundary term. 

This work has shown that & = V x p, which implies 
that the mean-field variation of p is given by p(r) = 
& X r in the continuum, up to an irrelevant constant of 
integration. If we choose p 9 = a x r 9 , then, on average, 
each p 9 is in the middle of the tile corresponding to g, as 
shown in Figure [3] The sums in (|22|) therefore scale with 
the number of terms, i.e, as Nrc ~ N, and Nec ~ V~N, 
respectively, and hence for large enough packings, the 
boundary term is negligible. In fact, with this choice of 



Polydisperse Packings 

For clarity we have considered only packings of identi- 
cal d-spheres, but many of our results extend to polydis- 
perse packings. Indeed, because it is topological in na- 
ture, the intrinsic formulation ([8]), (fT0|). and (fTT]) holds 
for packings of arbitrary convex grains, without modifi- 
cation. 

To extend the extrinsic formulation to polydisperse 
spheres, it is most convenient to use the radical Voronoi 
tesselation [20[ , which ensures that the edges of Voronoi 
cells pass through contacts. Discrete derivatives can be 
defined exactly as in this paper, and the definition of p is 
easily made, so that & = V x p. However, ip as defined in 
this paper does not describe identically torque-balanced 
configurations. The extension of ip to the polydisperse 
case will be discussed in future work. 
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Force Law 

For simplicity, the discussion above was limited to hard 
disks and spheres, interacting without cohesion. In fact, 
these assumptions are inessential. 

We considered hard disks and spheres, and framed this 
paper as a discussion of force configurations on a fixed 
geometry, to show the power of degree-of-freedom count- 
ing. However, the definitions of p and ip hold for disks 
and spheres of arbitrary softness. Provided contacts ap- 
pear at the midpoint of contact vectors, the definition of 
ip also holds. 

We also considered purely repulsive packings. In a re- 
pulsive packing, all real grains form closed loops. When 
grains admit cohesive forces, in addition to the loop- 
forming grains which bear the external load, there may 
be tree-like subsets of grains which extend into loops. 
Starting from the leaves of the tree, it is easy to see that 
each contact force on each of these grains must vanish. 
Therefore, provided the grains are treated as virtual, all 
of the results of this paper hold without modification. 
The physical novelty of cohesion is that the vanishing 
contact forces may be composed of an adhesive compo- 
nent, for example due to hydrodynamic forces, as well 
as an elastic component. This allows the grains to re- 
main in place under an infinitesmal disturbance, unlike 
the rattlers. 

More generally, this discussion shows that our results 
have little to do with the form of the force law and the 
repulsive constraints at the contacts, although the latter 
are an essential feature of granular physics [2l|. Rather, 
the key physical requirement is that grains interact lo- 
cally, so that a well-defined contact network exists. Our 
results then follow from the general form of Newton's 
laws, (JTJ), acting on this network. Our results may there- 
fore be useful in other problems involving local balance 



constraints, e.g. rigidity percolation [22| and metabolic 



networks [23 1 



CONCLUSION 



We have shown that Newton's laws ([1) can be inter- 
preted as differential relations in a discrete calculus. In 
an intrinsic topological formulation, these take the form 



d*/ = 0, d*(r x /) = 0. 



(26) 



These fields can be written in terms of gauge fields p and 

ip as 

f = d*p, rxf = d*(cp + rxp), (27) 
which are required to satisfy, for consistency, 

d*(<p + r x p) = r x d*p. (28) 



As discussed above, the stress tensor does not appear in 
this formulation, so the relation §4§ [ s absent. Equations 
(|26l) . ([27"!) and (|28| are exact, and hold for polydisperse, 
cohesive or non-cohesive, soft d-sphere packings. In this 
paper we considered d = 2 (disks) and d = 3 (spheres), 
but their extension to higher dimensions is straightfor- 
ward. 

To make contact with continuum equations, we also 
constructed an extrinsic formulation using the Delaunay 
triangulation and Voronoi tesselation. Using the triangu- 
lation, we constructed a discrete calculus in which New- 
ton's laws JT]) reproduce their continuum equivalent J3]) 
at the scale of a single grain, i.e., 



(V-<x) 9 = 0, 



= (a 9 



(29) 



In 2D, we introduced gauge fields p and ip, defined so 
that 



cr 9 = (V x p) 9 



p* = (V x 



(30) 



The same relations hold in 3D with p > p and ip — > 
To obtain physical force configurations, these must be 
subject to the geometry-dependent virtual contact con- 
straint ct c \ c£ vc = 0. Equations (f29|) . and ([30| are ex- 
act, and hold for monodisperse, cohesive or non-cohesive, 
packings of soft disks or spheres. 

Together, (30]) imply that j» = (VxVxfl», which 
is an exact, discrete representation of (4}. On a homoge- 
neous packing geometry, in the continuum limit, <x will 
satisfy <r = V x V x ip, with ip — > ip in 3D. From this 
expression we see that the pressure 



P S ±tr(*) 



iV 2 tr(t/;) - | V • (V • ^) 



in 2D 
in 3D. 



(31) 



We also proved that, in 2D, if macroscopic stresses are 
repulsive, then ip is a convex function. In general, these 
relations are insufficient to completely specify the con- 
tinuum limit of this problem, just as in elasticity a miss- 
ing equation is required that derives from Saint Venant's 
compatibility condition, via Hooke's law pij ]. 

The discrete representation of ip developed in this work 
allows insight into this missing stress-geometry equation 
d M, E3- Indeed, at the microscopic level, ip is only 
required to satisfy the virtual contact constraints, and 
the repulsive constraints. The former are present in all 
packing geometries but a very special one: the triangu- 
lar lattice in 2D. On this lattice, the grains are locally 
and globally as densely packed as possible. In a sense, 
the geometry is trivial. If we assume that the repulsive 
constraint requires only that macroscopic stresses are re- 
pulsive, then ip can be any convex function. 

In general, however, virtual contact constraints exist 
and couple the ip field throughout the packing. Their 
distribution is intimately related to the size and shape of 
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loops. For hyperstatic packings, Newton's laws are insuf- 
ficient to fully specify the stress state at the microscopic 
level, so we expect a family of solutions at the macro- 
scopic level as well. These could, in principle, depend on 
the microscopic force law and the packing history. In the 
simplest case, they would depend only on the distribution 
of virtual contacts. 

Only at isostaticity is the microscopic stress state fully 
specified by the geometry. In the special case of isotropic 
forcing, isostaticity is achieved in a noncohesive packing 
when the pressure P = 0. Using the microscopic ex- 
pression ([2]), we see that this implies J2 c grc I^ c |/at = 0, 
where is the normal component of the contact force 
at c. For non-cohesive grains, each /j^ > and hence 
all normal forces vanish. Assuming a Coulomb repulsive 
constraint of the form 

-\fr\<fh 
M 

for any /i < oo this implies that each tangential force 
vanishes as well, and hence the stress state is trivial. To 
understand the isostatic state for fj, < oo, forces need to 
be renormalized. 

However, for the ideal case fj, = oo, at P = only 
the normal components of all contact forces vanish. A 
nontrivial stress state is possible, consisting only of tan- 
gential forces. From (|31[) . we see that V 2 ip = in 2D 
and V 2 tr(i/j) = V • (V • tp) in 3D. In 2D, given boundary 
conditions on ip, this fully specifies the continuum limit. 
The stress state described by a = —Wip, with V 2 t/> = 
is equivalent to that of a linear elastic body undergoing 
volume-conserving deformation. 

To understand what happens at isostaticity when fi < 
oo, and when P > 0, requires more sophisticated anal- 
ysis. Since a noncohesive packing at isostaticity has no 
intrinsic force scale, by dimensional analysis, the missing 
equation must take the general form, in 2D, 



where T = \J P 2 — det(<x) is the maximal shear stress, 
and T is an unknown function of the geometry. In a 
forthcoming paper [26}, we use the tools of this work to 
derive the missing function T, using Edwards' statistical 
mechanics. 



APPENDIX A. TORQUE BALANCE IN 2D 

To verify that torque balance is satisfied identically 
when p = V x xp, we need to show that (V ■ p) 9 = 0. 
Labelling the contacts and triangles around a grain as in 




FIG. 4. Geometry used in torque balance computation. 



Figure |4j we note that s* is parallel to l c °^,t) ; so t\y & t 
tern 

= E 4 x ^ {-^r +f + f + +f°f°) 
= E ^ - A*r) 

t£Ts 
= 

on summation around the grain. 

APPENDIX B. DEGREES OF FREEDOM IN 3D 

Applying Eulcr's formula to the Voronoi polyhedron 
of grain g, we have Ny — + = 2. Each vertex 
is the confluence of three edges of the polyhedron, so 
that 3N V = 2Nj,. Summing over the entire packing, 
4AV - N B v - 3iV T + N BT + 2N C - N BC = %N and 
12AV — 3Nbv = 6-/Vt — 2Nbt- Since the entire packing 
is itself a convex polyhedron with trivalent vertices, we 
have Nbv-N B t + Nbc = 2 and 3N B v = 2N B t, so that 
4 Ay - 3 A T + 2 A c = 2 A + 2 and 2 Ay = A T . With these 
relations we can write M3 = 3 At — 3 Arc — 3 (Aye — 
Aye) — 3 (Ay — 1). This shows that the loop forces p are 
constrained by torque balance and the virtual contact 
constraints, Aye? of which are redundant. The gauge 
freedom has a dimension of Ay — 1 because adding a 
constant to B v does not change Ap. 

For V, we write M 3 = N c - 3Ay C + 6Ay G - 3(A T - 
A c ) - (A c - A) - (6A C - 6 - 12Ay + A). All but the 
first three terms correspond to gauge freedoms. First, we 
have gauge transformations of the form Ai/j c — r c ■ _B, 
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-8* x B. Here B is a con- 



which lead to (Ap)* c 

stant, but it may be derived from a fluctuating field on 
the triangles, viz. B = X)teT c ^ ■ This gives 3iV^ un- 
knowns constrained by 3Nq equations. Second, we can 
have Aip c = {i^l ■ (Ai/> 9+ - &ip a ~), with A^ 9 = r^B, 
where again B is a constant. In this case it can be de- 
rived from B = ^2 ceC g B c , which gives Nc unknowns 
constrained by N equations. Gauge transformations of 
this type leave invariant p. There remains a gauge sub- 
space of dimension 6Nc — 6 — 12Ny + N, the significance 
of which is unclear. 



APPENDIX C. CONVEXITY OF ip 

The force acting on a plane with unit normal n and 
area A is An ■ <r, with normal component An a n. This 
is always positive, and hence repulsive, if and only if <r 
is a positive-definite tensor. 

In 2D, <r = V x V x tjj can be inverted to write VVV' = 
— e ■ <r • e. Using a matrix identity, this can be written 



= tr(<x)5 - <r, 



(32) 



where 6 is the identity tensor. 

Writing H = VV^), convexity of ip is equivalent to 
positive definitcness of H. The latter is equivalent to 
the statement that H has positive eigenvalues. Writing 
Aj and Ui for the eigenvalues and eigenvectors of &, we 
have tr(o-) = Aj. From ([32]) . H ■ Uj — J2i^j ^i u j> so 
Uj is an eigenvector of H with eigenvalue Ylj-j-j Aj. This 
is positive if <r is positive definite. 
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